* Data analysis
	version 15
	capture ssc install coefplot
	
* Declare panel and ensure sort
	xtset ein year_soi
	sort ein year_soi
	
* Nix IVs for 2013-2019 for safety
	recode years_of_na_dummy fundr_ratio_dummy admin_ratio_dummy profit_ratio_dummy int_exp_ratio_dummy hhi_abs_dummy (nonmissing = .) if year_soi > 2012


* Fit models and store estimates; Also determine which cases the models run on
	* Full model for log_total_expenses
		xtreg log_total_expenses L.years_of_na_dummy L.fundr_ratio_dummy L.admin_ratio_dummy L.profit_ratio_dummy L.int_exp_ratio_dummy L.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 
		estimates store L1
		gen sample_L1 = e(sample)

		xtreg log_total_expenses L2.years_of_na_dummy L2.fundr_ratio_dummy L2.admin_ratio_dummy L2.profit_ratio_dummy L2.int_exp_ratio_dummy L2.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 	
		estimates store L2
		gen sample_L2 = e(sample)

		xtreg log_total_expenses L3.years_of_na_dummy L3.fundr_ratio_dummy L3.admin_ratio_dummy L3.profit_ratio_dummy L3.int_exp_ratio_dummy L3.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 	
		estimates store L3
		gen sample_L3 = e(sample)

		xtreg log_total_expenses L4.years_of_na_dummy L4.fundr_ratio_dummy L4.admin_ratio_dummy L4.profit_ratio_dummy L4.int_exp_ratio_dummy L4.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 	
		estimates store L4
		gen sample_L4 = e(sample)

		xtreg log_total_expenses L5.years_of_na_dummy L5.fundr_ratio_dummy L5.admin_ratio_dummy L5.profit_ratio_dummy L5.int_exp_ratio_dummy L5.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 	
		estimates store L5
		gen sample_L5 = e(sample)

		xtreg log_total_expenses L6.years_of_na_dummy L6.fundr_ratio_dummy L6.admin_ratio_dummy L6.profit_ratio_dummy L6.int_exp_ratio_dummy L6.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 	
		estimates store L6
		gen sample_L6 = e(sample)

		xtreg log_total_expenses L7.years_of_na_dummy L7.fundr_ratio_dummy L7.admin_ratio_dummy L7.profit_ratio_dummy L7.int_exp_ratio_dummy L7.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 	
		estimates store L7
		gen sample_L7 = e(sample)

		xtreg log_total_expenses L8.years_of_na_dummy L8.fundr_ratio_dummy L8.admin_ratio_dummy L8.profit_ratio_dummy L8.int_exp_ratio_dummy L8.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 	
		estimates store L8
		gen sample_L8 = e(sample)

		xtreg log_total_expenses L9.years_of_na_dummy L9.fundr_ratio_dummy L9.admin_ratio_dummy L9.profit_ratio_dummy L9.int_exp_ratio_dummy L9.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 	
		estimates store L9
		gen sample_L9 = e(sample)

		xtreg log_total_expenses L10.years_of_na_dummy L10.fundr_ratio_dummy L10.admin_ratio_dummy L10.profit_ratio_dummy L10.int_exp_ratio_dummy L10.hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 
		estimates store L10
		gen sample_L10 = e(sample)

* Use coefplot to make graphs
	graph drop _all
	* Years of net assets
	coefplot L1 L2 L3 L4 L5 L6 L7 L8 L9 L10, keep(L*.years_of_na_dummy) xline(0) nolabel legend(off) title(Years of Net Assets Coef.) ytitle(Lag (in years)) ylabel(1 2 3 4 5 6 7 8 9 10) msymbol(O) scheme(sj) nooffsets note(Note: Coef. indicates whether years of net assets is above the sector-year median) name(yona)
* Fundraising ratio
	coefplot L1 L2 L3 L4 L5 L6 L7 L8 L9 L10, keep(L*.fundr_ratio_dummy) xline(0) nolabel legend(off) title(Fundraising Expense Ratio Coef.) ytitle(Lag (in years)) ylabel(1 2 3 4 5 6 7 8 9 10) msymbol(O) scheme(sj) nooffsets note(Note: Coef. indicates whether fundraising expense ratio is above the sector-year median) name(fundr)
	* Admin ratio
	coefplot L1 L2 L3 L4 L5 L6 L7 L8 L9 L10, keep(L*.admin_ratio_dummy) xline(0) nolabel legend(off) title(Administrative Expense Ratio Coef.) ytitle(Lag (in years)) ylabel(1 2 3 4 5 6 7 8 9 10) msymbol(O) scheme(sj) nooffsets note(Note: Coef. indicates whether administrative expense ratio is above the sector-year median) name(admin)
	* Profit ratio
	coefplot L1 L2 L3 L4 L5 L6 L7 L8 L9 L10, keep(L*.profit_ratio_dummy) xline(0) nolabel legend(off) title(Profit Ratio Coef.) ytitle(Lag (in years)) ylabel(1 2 3 4 5 6 7 8 9 10) msymbol(O) scheme(sj) nooffsets note(Note: Coef. indicates whether profit ratio is above the sector-year median) name(profit)
	* Interest expense ratio
	coefplot L1 L2 L3 L4 L5 L6 L7 L8 L9 L10, keep(L*.int_exp_ratio_dummy) xline(0) nolabel legend(off) title(Interest Expense Ratio Coef.) ytitle(Lag (in years)) ylabel(1 2 3 4 5 6 7 8 9 10) msymbol(O) scheme(sj) nooffsets note(Note: Coef. indicates whether interest expenses ratio is above the sector-year median) name(interest)
	* HHI
	coefplot L1 L2 L3 L4 L5 L6 L7 L8 L9 L10, keep(L*.hhi_abs_dummy) xline(0) nolabel legend(off) title(Herfindahl-Hirschman Index (HHI) Coef.) ytitle(Lag (in years)) ylabel(1 2 3 4 5 6 7 8 9 10) msymbol(O) scheme(sj) nooffsets note(Note: Coef. indicates whether HHI is above the sector-year median) name(hhi)

* Obtain cumulative effect estimates in one omnibus model
	xtreg log_total_expenses L(1/10).years_of_na_dummy L(1/10).fundr_ratio_dummy L(1/10).admin_ratio_dummy L(1/10).profit_ratio_dummy L(1/10).int_exp_ratio_dummy L(1/10).hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, fe vce(robust) 
	gen sample_Omni = e(sample)

	* Postestimation: Get cumulative marginal effects for each DV over the lags
		* Years of net assets:
		local yona_effect = _b[L.years_of_na_dummy] + _b[L2.years_of_na_dummy] + _b[L3.years_of_na_dummy] + _b[L4.years_of_na_dummy] + _b[L5.years_of_na_dummy] + _b[L6.years_of_na_dummy] + _b[L7.years_of_na_dummy] + _b[L8.years_of_na_dummy] + _b[L9.years_of_na_dummy] + _b[L10.years_of_na_dummy]
		display `yona_effect'
		* Fundraising expense ratio: 
		local fundr_effect = _b[L.fundr_ratio_dummy] + _b[L2.fundr_ratio_dummy] + _b[L3.fundr_ratio_dummy] + _b[L4.fundr_ratio_dummy] + _b[L5.fundr_ratio_dummy] + _b[L6.fundr_ratio_dummy] + _b[L7.fundr_ratio_dummy] + _b[L8.fundr_ratio_dummy] + _b[L9.fundr_ratio_dummy] + _b[L10.fundr_ratio_dummy]
		display `fundr_effect'
		* Administrative ratio:
		local admin_effect = _b[L.admin_ratio_dummy] + _b[L2.admin_ratio_dummy] + _b[L3.admin_ratio_dummy] + _b[L4.admin_ratio_dummy] + _b[L5.admin_ratio_dummy] + _b[L6.admin_ratio_dummy] + _b[L7.admin_ratio_dummy] + _b[L8.admin_ratio_dummy] + _b[L9.admin_ratio_dummy] + _b[L10.admin_ratio_dummy]
		display `admin_effect'
		* Profit ratio:
		local profit_effect = _b[L.profit_ratio_dummy] + _b[L2.profit_ratio_dummy] + _b[L3.profit_ratio_dummy] + _b[L4.profit_ratio_dummy] + _b[L5.profit_ratio_dummy] + _b[L6.profit_ratio_dummy] + _b[L7.profit_ratio_dummy] + _b[L8.profit_ratio_dummy] + _b[L9.profit_ratio_dummy] + _b[L10.profit_ratio_dummy]
		display `profit_effect'
		* Interest expense ratio:
		local interest_effect = _b[L.int_exp_ratio_dummy] + _b[L2.int_exp_ratio_dummy] + _b[L3.int_exp_ratio_dummy] + _b[L4.int_exp_ratio_dummy] + _b[L5.int_exp_ratio_dummy] + _b[L6.int_exp_ratio_dummy] + _b[L7.int_exp_ratio_dummy] + _b[L8.int_exp_ratio_dummy] + _b[L9.int_exp_ratio_dummy] + _b[L10.int_exp_ratio_dummy]
		display `interest_effect'
		* HHI:
		local hhi_effect = _b[L.hhi_abs_dummy] + _b[L2.hhi_abs_dummy] + _b[L3.hhi_abs_dummy] + _b[L4.hhi_abs_dummy] + _b[L5.hhi_abs_dummy] + _b[L6.hhi_abs_dummy] + _b[L7.hhi_abs_dummy] + _b[L8.hhi_abs_dummy] + _b[L9.hhi_abs_dummy] + _b[L10.hhi_abs_dummy]
		display `hhi_effect'
		display "Total of cumulative effects: " `yona_effect' + `fundr_effect' + `admin_effect' + `profit_effect' + `interest_effect' + `hhi_effect'

* Obtain summary statistics 
	* Determine which cases the models run on to obtain relevant summary statistics
		egen newvar = rowtotal(sample_L1 sample_L2 sample_L3 sample_L4 sample_L5 sample_L6 sample_L7 sample_L8 sample_L9 sample_L10 sample_Omni)
		gen case_used = 0
		recode case_used (0=1) if newvar >0 & newvar < .
		drop newvar
		label var case_used "Indicator that the case was used in at least one model estimation"
		tab case_used

* Obtain counts
	count
	count if year_soi == .
	count if year_soi != . & out_nccs == "IN" & donative != 1
	count if year_soi != . & out_nccs == "IN" & donative == 1
	count if year_soi != . & out_nccs == "IN" & donative == 1 & case_used == 1
	count if year_soi != . & out_nccs == "IN" & donative == 1 & case_used != 1
	
	bysort year_soi: count
	bysort year_soi: count if case_used==1
		
* Obtain summary statistics restricted by case_used
	tabstat pct_donative total_expenses admin_ratio fundr_ratio hhi_abs years_of_na profit_ratio int_exp_ratio if case_used==1, columns(statistics) statistics(mean sd median count min max) varwidth(20) format 
	tab ntmaj12 if case_used==1

* Check
	bysort year_soi: tabstat pct_donative total_expenses admin_ratio fundr_ratio hhi_abs int_exp_ratio profit_ratio years_of_na if case_used==1, columns(statistics) statistics(mean sd median count min max) varwidth(20) format
	* The IRS SOI extracts do not include data on administrative and fundraising expenses or government grants after 2012, so data on the administrative expense ratio, fundraising expense ratio, and the HHI are missing for those years. Cases for 2013-2019 are included in lagged models when the independent variables are available in years prior to 2013. 

	sort ein year_soi

*Address Endogeneity 
	reg log_total_expenses D.L(1/10).years_of_na_dummy D.L(1/10).fundr_ratio_dummy D.L(1/10).admin_ratio_dummy D.L(1/10).profit_ratio_dummy D.L(1/10).int_exp_ratio_dummy D.L(1/10).hhi_abs_dummy i.year_soi if out_nccs=="IN" & donative==1, vce(robust) 
	*Lagged first difference model to address unobserved heterogeneity and reverse causality

* restore

